# This code is hosted on http://code.google.com/p/lenthorp/
# Freely available for use in applications, but should NOT be modified
# Email all comments to lenthorpresearch@gmail.com

import numpy as np
import csv

class Node:

	def __init__(self,centre,type):
		self.centre = centre
		self.type   = type
	def __init__(self,centre):
		self.centre = centre
		self.type   = 0
		
		
class TimeGrid:

	def __init__(self,endTime,startTime,nPts):
		self.nPts      = nPts
		self.endTime   = endTime
		self.startTime = startTime
		self.dt        = np.zeros(nPts - 1)
		self.timePts   = np.zeros(nPts)
		self.buildGrid()
		
	def buildGrid(self):
		dtt = (self.endTime - self.startTime)/(self.nPts - 1)
		self.timePts = np.array([self.startTime + i * dtt for i in range(self.nPts)])
		self.dt      = np.diff(self.timePts)
		
	def timeIndexAfterEqualT(self,time):
		i = -1
		for i in self.timePts:
			if self.timePts >= time : 
				return i
		return i

class RBFPricingEngine:

	def __init__(self,nodes,eqnL,boundaryConditions,timeGrid):
		self.nodes = nodes
		self.eqnL  = eqnL
		self.numNodes = len(nodes)
		self.timeGrid = timeGrid
		self.numInteriorNodes = 0
		self.solution = np.zeros(self.numNodes)
		self.boundaryCondtions = boundaryConditions
		self.alphas = []
		
		for n in self.nodes:
			if n.type == 0 :
				self.numInteriorNodes+=1
				
		self.interiorNodeIndicies = []
		for i in range(len(nodes)):
			if nodes[i].type == 0:
				self.interiorNodeIndicies.append(i)
				
		self.boundaryNodeIndicies = []
		for i in range(len(nodes)):
			if nodes[i].type != 0:
				self.boundaryNodeIndicies.append(i)
		
	def calculateAtNode(self,rbfFunc,node)		:
		return np.sum(rbfFunc(node,self.nodes[i])*self.alphas[i] for i in range(len(self.nodes)))
	
	def solveConstCoef(self,fromTime,toTime,initialSolution,rbfZeroDeriv):
		##build equation to solve
		n = len(self.nodes)
		D = np.zeros([n,n])
		A = np.zeros([n,n])
		self.solution = initialSolution

		##interior nodes points
		for i in self.interiorNodeIndicies:
			nd = self.nodes[i]
			A[i] = np.array([rbfZeroDeriv(nd,xc) for xc in self.nodes])
			for x in self.eqnL:
				func = x[1]
				coef = x[0]
				D[i] += np.array([coef * func(nd,xc) for xc in self.nodes])
			
				
		for i in self.boundaryNodeIndicies:
			nd = self.nodes[i]
			bc = self.boundaryCondtions[nd.type]
			func = bc[0]
			D[i] = np.array([func(nd,xc) for xc in self.nodes])
			self.solution[i] = -bc[1]
			
		
		indexFrom = np.nonzero(self.timeGrid.timePts>= fromTime)[0][0]
		indexTo   = np.nonzero(self.timeGrid.timePts<= toTime)[0][0]
		time      = fromTime
		nk        = indexFrom - indexTo 
		
		for k in range(nk):
			dt = self.timeGrid.dt[k]
			
			self.solution[self.interiorNodeIndicies]/=dt
			self.alphas  = np.linalg.solve(D-A/dt,-self.solution)

			##update solution
			self.solution = np.array([np.sum(rbfZeroDeriv(x,self.nodes[j])*self.alphas[j]
											 for j in range(n)) for x in self.nodes])
			for i in self.boundaryNodeIndicies:
				self.solution[i] = -self.boundaryCondtions[self.nodes[i].type][1]

		self.solution = np.array([np.sum(rbfZeroDeriv(x,self.nodes[j])*self.alphas[j]
											 for j in range(n)) for x in self.nodes])
		

			
			
			
		
		

			
		








		